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UNSTEADY  THERMAL  BLOOMING  OF  INTENSE  LASER  BEAMS 


INTRODUCTION 

When  laser  light  passes  through  an  absorbing  fluid,  the  beam  intensity  is  reduced  by 
molecular  absorption.  If  the  fluid  contains  inhomogeneities  or  impurities  (for  example, 
particulate  matter  or  turbulent  eddies  in  the  case  of  a  gas)  other  scattering  and  absorption 
mechanisms  can  contribute  to  intensity  reduction.  This  report  is  confined  to  a  study  of 
laser  beam  propagation  through  an  absorber  for  which  it  is  appropriate  to  describe  the 
interaction  of  the  light  and  the  fluid  by  means  of  an  absorption  constant  per  unit  length 
of  absorbing  path.  At  low  powers  a  single  absorption  constant  a  suffices  to  describe  the 
beam  diminution  even  in  the  presence  of  scattering  centers.  At  high  powers,  however,  the 
absorbed  energy  causes  significant  enough  alteration  of  the  ambient  fluid  density  that  the 
light  beam  itself  is  perturbed,  not  simply  diminished.  This  phenonenon  of  self-distortion 
of  the  coherent  phase  front  due  to  absorber  heating  is  called  thermal  blooming.  If  scattering 
is  occurring  as  well,  then  it  is  important  to  separate  the  purely  “elastic”  losses  of  the  beam 
via  scattering  from  the  “unelastic”  losses  via  absorption  so  that  the  absorption  coefficient 
is  written 


«  "  <W  +  ««*• 

Thermal  blooming  has  been  studied  under  a  variety  of  limiting  conditions.  If  a  wind 
is  present  in  the  absorbing  gas,  or  if  the  beam  itself  is  moving  through  the  fluid,  then  the 
heated  volume  of  the  absorber  is  removed  from  the  beam  path  after  residing  a  finite  time 
(about  equal  to  beam  size/fluid  speed  relative  to  the  beam)  and  a  steady  state  can  evolve. 
Time  no  longer  explicitly  enters  the  problem  and  a  three-space-dimensional  analysis  is 
appropriate.  The  resultant  nonlinear  problem  has  been  studied  exhaustively  at  NRL  and  a 
number  of  other  laboratories,  [1]  Another  limit  of  interest  is  the  regime  of  times  which  is 
much  shorter  than  the  wind  transit  time  across  the  beam,  Then,  if  the  initial  beam  is  dis¬ 
tributed  symmetrically  about  the  propagation  axis,  the  thermal  blooming  evolves  sym¬ 
metrically  about  the  propagation  axis,  the  thermal  blooming  evolves  symmetrically  and, 
again,  just  three  coordinates  are  needed:  the  radius,  the  range,  and  time.  This  regime  also 
has  been  studied.  [1] 

There  are,  however,  other  space-time  regimes  of  interest  and  importance  in  naval 
applications  of  high-energy  lasers,  such  as  arbitrary  (noncircularly  symmetric)  aperture 
distributions  in  the  short-time  regime,  arbitrary  distributions  in  the  regime  before  steady 
state  is  attained,  propagation  through  wind-null  points  in  the  beam  path,  propagation  of 
time-varying  aperture  distributions,  and  multiple  pulse  train  propagation,  These  require  a 
four-dimensional  analysis. 
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This  report  details  the  development  of  a  four-dimensional  computer  program  which  is 
being  used  to  study  the  preceding  problems.  The  physics  of  the  light  beam,  the  fluid,  and 
their  interaction  are  discussed.  A  section  is  devoted  to  some  numerical  details,  a  comparison 
with  some  experimental  data  is  given,  and  a  flowchart  is  provided  as  an  appendix. 


PHYSICS  OF  THE  LASER-FLUID  INTERACTION 

When  the  laser  is  first  turned  on ,  a  transient  state  evolves  during  which  the  tempera¬ 
ture  of  the  fluid  increases  due  to  local  heating  at  the  rate  oi^Hx,  y,  z,  0  ergs  per  unit 
volume  per  unit  time,  where  /  is  the  laser  beam  intensity.  A  pressure  wave  propagates  away 
from  the  beam  at  the  speed  of  sound.  Directed  velocity  of  the  fluid  away  from  the  beam 
region  leads  to  density  decrease  there.  This  study  is  confined  to  times  after  these  transients 
have  died  out,  which  times  are  on  the  order  of  the  beam  size/sound  speed.  In  this  limit, 
as  will  be  shown  below,  the  density  changes  are  proportional  to  the  heat  absorbed  by  a  fluid 
element  during  its  residence  in  the  beam. 

The  fluid  equations  to  be  solved  are  conservation  of  mass 


(1) 


conservation  of  momentum 


pgf  +  »  0,  (2) 

and  conservation  of  energy 

9 |f  +  pV«v  -  Q  ,  (3) 


where 


J?  =  ! 

Dt  dt 


+ 


vV. 


In  the  equations  p ,  v,  and  p  are  the  density,  velocity,  and  pressure,  e  is  the  internal  energy 
per  unit  volume  e  ■  CyT ,  where  Cv  is  the  specific  heat  at  constant  volume  per  unit  volume, 
T  is  the  temperature,  and  Q  is  the  rate  of  heat  deposition  by  the  laBer  beam  per  unit  volume, 
so  that  Q  *  otabll  I, 

In  the  time  regime  in  the  preceding  discussion  pressure  is  considered  constant.  In 
addition,  an  ideal  gas  equation  of  state  is  assumed.  Then  Eqs.  (1),  (2),  and  (3)  become 

Dt  ”ix  (I'-1)CW’  <4> 
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where  T  =  Cp/Cv  and  Cp  is  the  specific  heat  at  constant  pressure  per  unit  volume,  and  pu  is 
the  ambient  pressure. 

In  all  cases  of  interest  Ap  0  (a)  and  linearization  of  the  hydrodynamics  is  performed. 
This  assumption  must  be  checked  in  some  cases  (for  example,  long  pulses  of  high  energy 
in  a  confined  absorption  cell  can  lead  to  significant  temperature  changes  compared  with 
ambient  values  [2] ).  For  naval  applications,  however,  the  assumption  is  satisfied  very  well. 
Thus, 


p  -  p0  +  pj.pj  «pQ 


and  Eq.  (4)  becomes 


DPl 

Dt 


IEjlU 

r  2 


<w  - 


where  Cg  is  the  sound  speed, s/Pp0/p0.  The  solution  to  Eq.  (6)  is 


(5) 


Pi  (*i  2,  0  ■  otabl  j  /[»  -  v(f  -  t  ]  dt (6) 

CS  J0 

This  result  reduces  to  some  of  the  limiting  cases  which  were  the  subject  of  earlier  studies,  as 
will  now  be  shown. 

For  v  -  0,  and  short  times  (longer,  however,  than  the  beam  size/sound  speed)  Eq.  (6) 
becomes 


p^r.z,  t)  *»  -(r-l)aabs/(r,  z,  t)t, 


(7) 


which  is  the  density  change  appropriate  to  long  pulses  for  which  negligible  degradation  has 
occurred  .during  the  transient  sound  wave  period  [3]  (r  -*■  |r|  for  circularly  symmetric  initial 
conditions). 

For  v  0  anywhere  between  aperture  and  focal  point,  and  for  (  -*  Eq.  (6)  becomes, 
with  r'  ■  r  -  v(f  -  t  ’) 


Pi  (r,  z) 


-  i)«ab,  rr 

J-c 


Icw(r' ,z)dt' , 


(8) 


which  is  the  form  of  the  density  changes  in  the  CW  propagation  studies  (4  j . 
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in  the  current  study,  Eq.  (6)  is  used  without  further  simplification,  In  its  present  form, 
it  is  inapplicable  to  cases  where  |v|  approaches  the  sound  speed,  so  that  the  pressure  waves 
are  always  in  the  vicinity  of  the  beam.  Modifications  to  account  for  high  velocities  are  easily 
accomplished  [5]  and  could  be  included  for  future  studies. 

The  laser  beam  is  described  by  the  scalar  wave  equation  in  the  paraxial  approximation. 
The  derivation  of  this  equation  is  detailed  elsewhere  [4]  and  Is 


+  ^  +  lii  + 

bx2  by2 


k2(ri2  -l)i p 


0, 


(9) 


where  is  the  complex  scalar  ws1  e  amplitude  (electrL  iield  E  *  exp [i(kz  -  ltd)) ),  k  is  the 
radian  wavenumber  2ir /X,  and  n(x,  y,  z,  t)  is  the  index  of  refraction,  with  its  ambient  value 
being  taken  to  be  n  ”  1.  Equation  19)  holds  for  slowly  varying  amplitude  in  the  z  direction, 

I  kb\p/b z  |  »  I  b2\p/b z2  |.  The  intensity  is  C/8ir  !  \p  |2e*°*.  Equation  (t>)  is  coupled  to  Eq. 
(6)  via  the  Lorentz-Lorentz  law, 


n2-l 
n2  +  2 


Np, 


(10) 


where  N  is  the  molecular  refractivity 


N 


2 

3  P0 


bn 

bT 


(ID 


for  a  perfect  gas,  ambient  index  equal  to  one,  and  small  density  changes  p, , 

Note  that  the  time  dependence  in  Eq.  (9)  is  Implicitly  contained  In  the  change  in  index 
n2  -  1.  The  high-frequency  electric  field  oscillations  have  been  removed,  and  the  speed  of 
light  has  been  taken  to  be  infinite  (neglecting  retardation).  Thus  the  light  beam  is  altered  on 
a  time  scale  determined  by  the  fluid  motion. 


COMPUTATIONAL  DETAILS 

Equations  (6)  (9),  and  (10)  are  solved  in  the  following  sequence  (as  given  by  the 
flowchart  in  the  Appendix). 

The  amplitude  appropriate  for  the  apertu  o  distribution  of  interest  is  given  at  z  »  0  for 
all  t.  Using  these  values  the  density  change  b>  ween  z  n  0  and  z  -  Az  is  predicted  via  Eq.  (6) 
and  „he  index  change  via  Eq.  (10),  The  wave  amplitude  is  propagated  to  z  -  A z,  all  t ,  via  Eq. 
(9),  and  this  new  intensity  pattern  is  used  to  recalculate  the  index  field  to  correct  the 
amplitudes.  This  completes  the  basic  iteration  module  from  which  a  solution  is  constructed 
by  marching  from  aperture  to  focus  or  range  of  interest. 
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The  time  integration  in  Eq.  (6)  must  be  treated  with  care.  Since  the  integrand  depends 
on  the  upper  limit  it  appears  that  the  integral  must  be  recomputed  completely  at  each  time 
step.  This  can  be  avoided  as  follows:  The  density  at  t  +  A t  is,  by  definition 


p(r,  z,  t  +  At)  -  -  - 


c  2 


I  +  A  / 

7[r -  v(f  +  At  -  t'),  z,  t']  dt\  (12) 


and  we  wish  to  relate  it  to  p(r,  y,  f).  We  have 

(r-i) 


p(r,  z,  t  +  At)  -  - 


r  ^ 


aot»s  r 
Jo 


/(r-v(f  +  At  -  f'),z,  f  ’]  dt ' 


(r-D«atw 


I. 


t  +  Af 

/[r-v(t  +  Af-f’),z,  t']  dt' 


p(t  -  vAf,  z,  t)  - 


(r-l)OUa 


c  2 

‘-'s 


s. 


t  +  A< 


/|r  -  v(t  +  At  -  t  '),2,  t  ’]  dt ‘ 


The  integration  can  be  done  by  trapezoidal  rule  to  give 


<r-D«abs 


Af 


p(r,  z,  t  +  Af)  ■  p(r  -  vAf,  z,  t) - — ~  [/(r,z,  f  +  Af)  +  /(r  -  vAf,  z,  f))  ^  .  (13) 


We  do  not  know  /(r,  z,  f  +  At)  at  this  point  in  the  calculation,  since  we  need  p(r,  z,  f  +  Af) 
to  get  it.  Thus  we  estimate  /( r,  z,  f  +  At)  from  the  time  derivative 

I(t,  z,  f+Af)  /(r,  z,  t)  +  At  |“r  (r,  z,  f)  + - 

We  replace  the  time  derivative  by  a  central  difference 


9 1  I(t ,  z,  f+Af)  -  /(r,  z^f  -  At)  . 

at  . . 2Af 
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thus, 


J(r,  z.t+At)  =*  2/(r,  2,  I)  -  7(r,  2,  t -- At). 


Finally, 


(r-i)aD,„ 

p(r,  2, 1  +  At)  -  p(r  -  vAt,  2,  f) - 5-^  [2/(r,  2, 1) 

cs2 


At 


-  /( r,  2,  t  -  At)  +  I(r  -  vAt,  2,  t)] 


(14) 


FREE  CONVECTION 

For  beams  propagating  in  the  atmosphere,  the  buoyancy  induced  by  laser  heating  can 
be  neglected  [3] ,  In  laboratory  simulations,  however,  free  convection  can  produce  signifi¬ 
cant  flow  of  absorber  in  a  time  comparable  to  the  transient  heating  period  of  thermal 
blooming.  Buoyancy  should  not  be  ignored  in  the  experiments  UBed  to  check  the  computer 
code  and  in  the  following  discussion.  Thus  a  calculation  of  the  induced  vertical  velocity  was 
incorporated  in  the  code  in  order  to  achieve  a  satisfactory  quantitative  confirmation  of  the 
computer  calculation. 

Therefore,  in  addition  to  Eq.  (5),  we  must  also  solve,  neglecting  viscosity, 

/duB  \ 

pO\~dT  +  v'AvbJ  “  pl*<  <1B) 


where  g  is  a  negative  number,  representing  the  acceleration  of  gravity.  The  solution  to  Eq. 
(16)  is 


vB(x,y,  t)  -  -A 


r 


Pi[*  -  Vv(f  -  t ' ),  y  -  Mf  -  t'),  f']  df 


(16) 


And  due  to  the  presence  of  a  vertical  component  of  velocity,  Eq.  (6)  becomes 

Pi(*.  y.  *)  -  f  V*  -  t'),  t']  df. 

CS  Jo 


(17) 
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Combining  Eqs,  (15)  and  (16)  and  integrating  by  parts  gives 

"A  m^T  [  *‘nx-»x(t-t'),y-vB(t-f),f]  df.  (18) 

PoCs  Jo 

The  computer  code  must  solve  Eqs.  (17)  and  (18)  simultaneously.  This  is  accomplished  by 
alternating  the  calculation  of  p j  and  vB  and  iterating  until  each  converges.  In  practice,  it 
was  found  that  one  iteration  was  sufficient.  The  algorithm  for  calculating  vB  is  similar  to 
the  method  used  for  p1  and  given  in  Eq.  (14).  Define 

f  (19) 

Pocs2 


then  the  algorithm  is 

ufl(*,  y,  t  +  Af) 


8 px(0 

Po' 


(t+At)  +  J(x  -  vx£±t,  y  -  oB  At,  t) 


+  M  W  +  t2/(je,  y,  t)  -  /(*,  y,  t  -  At)) 

PoV  2  1 

+  tf(x  -  vxAt,  y  -  vB&t,  t)}  , 

and  Eq.  (14)  becomes 

p^x.yt  +  At)  -  p(x-vxAt,y~vBAt,t) 

■  (-“t-  t  A0 

cs 

+  /(x  -  u^At,  y  -  UgAf,  t)] . 


(20) 


(21) 


These  algorithms  require  the  storage  of  two  earlier  values  of  the  intensity  field,  /,  one 
earlier  value  of  the  integral  J  and  the  current  values  of  uB  and  p. 

Earlier  in  this  section  it  was  stated  that  the  equations  for  the  conservation  of  mass  and 
energy  would  be  solved  in  the  isobaric  limit,  This  will  now  be  justified  by  considering  the 
time  scales  of  importance  in  the  experiment  which  was  used  to  provide  validation  of  the 
computer  code  [6] .  A  pivoting  gas  absorption  cell  provided  a  null  convection  point  re¬ 
quiring  a  full  four-dimensional  analysis  (Pig.  1).  The  appropriate  time  interval  for  use  in  the 
computations  follows  from  these  considerations. 
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In  setting  up  the  requisite  hydrodynamic  equations  to  solve  the  propagation  problem, 
one  can  take  advantage  of  the  large  difference  between  time  scales  characterizing  pressure- 
wave  transients  and  convective  effects,  The  former  travel  at  the  speed  of  sound  are 
therefore  lost  from  the  beam  region  in  a  time  like  a{z)/cH,  which  is  about  10  ps  for  these 
experiments.  The  origins  of  convection  effects  are  twofold:  (1)  the  relative  velocity  between 
gas  and  laser  beam  due  to  pivoting  of  the  cell  and  (2)  the  induced  buoyancy  due  to  the 
heating  of  the  gas,  To  compare  the  time  scale,  due  to  cell  pivoting  with  the  pressure 
equilibration  time  calculated,  we  chose  the  maximum  velocity  induced  by  pivoting  in  all 
of  experiments  used  for  comparison.  This  occurred  when  the  cell  was  pivoted  about  the 
point  z  =  -  10  cm  with  a  slew  rate  of  0,2  rad/s,  giving  a  characteristic  time  (at  the  beam 
waist)  of  3.63  ms,  All  other  velocities  at  different  slew  rates,  ranges,  and  pivot  points  are 
slower.  Thus  the  forced  convection  times  are  at  least  three  orders  of  magnitude  slower  than 
the  time  characteristic  of  the  pressure  effects.  The  buoyancy-induced  velocities  were  esti¬ 
mated  and  found  to  be  greater  than  0,04  s,  and  therefore  comparable  to  thermal-blooming 
decay  times, 

The  time  scale  of  the  various  hydrodynamic  phenomena  suggests  that  we  can  ignore 
sound  effects  and  solve  merely  the  simpler  problem.  Beam  degradation  is  taken  to  be  and 
was  found  to  be  negligible  during  sound  transit  times.  In  addition,  the  effects  of  the  pressure 
waves  striking  the  wall  of  the  gas  cell  and  returning  to  the  region  of  the  beam  are  ignored 
and  were  suppressed  in  the  experiment  by  suitable  choice  of  all  walls, 

Great  care  was  taken  to  model  the  experiment  as  closely  as  possible.  The  beam  was 
propagated  without  absorption  from  z  a  0  to  z  80  cm,  At  2  *  80  crn  the  absorption 
coefficient  assumed  the  value  0.42  m"1 .  The  absorption  was  droppod  to  zero  at  103  cm 
farther  on,  and  the  final  values  were  calculated  10  cm  beyond  this  point.  The  nominal  beam 
sizes  of  0.32  cm  and  0,08  cm  at  entrance  and  exit  could  not  be  achievod  with  diffraction- 
limited  unperturbed  propagation,  Hence,  wavelength  scaling  was  used  to  produce  a  focal 
spot  1.43  times  the  diffraction  limit  to  fit  the  beam-size  constraint.  A  computational 
‘‘detector*’  was  used  by  counting  the  flux  at  mesh  points  which  fell  inside  the  circle  with 
radius  equal  to  the  radius  of  the  experimental  detector.  Thus  the  code  could  give  in  prin¬ 
ciple  an  absolute  comparison  with  the  experiment,  For  the  purposes  of  this  work,  however, 
all  flux  values  were  divided  by  the  initial  value  to  produce  data  similar  to  that  taken  ex¬ 
perimentally,  The  computational  detector  was  made  to  sense  the  peak  intensity  point  and 
followed  it  in  the  horizontal  cell  model,  The  detector  remained  on-axis, 

The  necessity  to  introduce  natural  convection  into  the  code  became  apparent  when 
the  horizontal  experiment  was  modeled  without  natural  convection.  The  data  indicated  too 
much  heating  gas  in  the  computational  model  and  suggested  the  presence  of  another  cooling 
mechanism  beside  the  forced  convection.  The  hydrodynamics  were  then  altered  to  produce 
the  algorithm  in  the  preceding  discussion. 


SOURCES  OF  ERROR 

The  following  approximations  have  contributed  to  the  inuccuracy  of  the  code  results: 

•  Wavelength  scaling  to  account  for  nondiffraction-limited  behavior.  This  gives  results 
that  are  most  likely  slightly  high, 
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•  Coarse  discretization  in  a,  y,  and  2  to  avoid  long  running  times.  This  can  give 
answers  either  too  high  to  too  low. 

•  Coarse  discretization  in  t  coupled  with  crude  (trapezoidal  rule)  time  integration.  The 
direction  of  this  error  is  unknown. 

•  Use  of  only  one  iteration  in  solving  Eqs.  (14)  and  (15).  This  leads  to  an  underestimate 
of  va  and  thus  lower  results  than  if  enough  iterations  for  convergence  were  used. 

•  The  approximation  of  Eqs.  (17)  and  (18).  This  is  a  predictor  without  a  corrector, 
and  since  'dl/bt  <  0,  it  gives  rise  to  a  lower  value  of  I.  Lower  values  of  I  will  give 
smaller  density  changes  and  smaller  buoyancy  velocities  which  tend  to  cancel  to 
order  At. 

If  one  assumes  all  of  these  effects  to  be  random  and  at  most  10%,  one  can  account  for 
an  rms  error  of  about  *  5%,  and  the  computer-generated  points  should  be  taken  to  fall 
within  these  ei.or  bounds. 


COMPARISON  WITH  EXPERIMENT 

In  the  experiment  [6] ,  a  C02  laser  beam  of  nominal  10-W  power  was  passed  through 
an  absorbing  cell  which  had  the  capability  of  being  pivoted.  The  pivot  point  was  variable, 

If  2g  is  the  position  of  the  pivot  point,  then  the  forced  convection  at  any  range  z  is  given  by 


(22) 


where  SI  is  the  pivot  rate.  If  the  pivot  point  is  inside  the  absorbing  cell,  there  is  no  forced 
convection  at  thatpoint,  and  buoyance  develops  as  the  sole  heat  removal  agent.  The  experi¬ 
ment  provided  a  good  test  of  the  computer  calculation. 

Figures  2  through  7  are  plots  of  the  relative  intensity  as  a  function  of  time  for  pivot 
points  varying  from  2  =  -  45  cm  to  zs  =  90  cm,  with  2  =  0  taken  at  the  entrance  of  the 
1 -meter  cell.  The  relative  intensity  is  defined  as 


^meas^^v 


o-a  L 


(23) 


where  /  is  the  reference  level  obtained  with  the  cell  evacuated,  a  is  the  absorption  coeffi¬ 
cient  of  the  gas  absorber,  and  L  is  the  length  of  the  cell,  The  solid  points  are  experimental, 
and  the  open  points  are  obtained  from  the  four-dimensional  code.  The  agreement  is  quite 
good  in  all  cases.  Improvement  could  be  achieved  by  finer  sampling  of  the  calculation  in 
range,  especially  in  the  vicinity  of  a  pivot  point  within  the  cell,  Since  the  code  calculations 
are  already  quite  time  consuming,  further  refinements  were  not  warranted  in  view  of  the 
satisfactory  results  obtained.  This  agreement  is  close  enough  for  the  purposes  for  which  the 
code  was  developed, 
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CONCLUSIONS 

A  four-dimensional  computer  code  has  been  written  which  solves  for  the  propagation 
of  an  intense  laser  beam  passing  through  an  absorbing  fluid.  Included  are  the  calculations 
of  complex  hydrodynamic  flows  induced  by  the  beam  itself,  The  code  was  confirmed  by 
comparing  its  results  with  those  of  an  experiment.  This  comparison  required  (for  agreement) 
that  the  code  include  ail  features  of  the  interaction  of  beam  and  absorber  on  each  other 
in  a  self-consistent  way. 
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